      subroutine parrefl(ncellsp, ijkctmp, iphead, npart, itdim, iwid, jwid, &
         kwid, x, y, z, wate, bxv, byv, bzv, px, py, pz, up, vp, wp, pxi, peta&
         , pzta) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncellsp 
      integer , intent(in) :: npart 
      integer  :: itdim 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ijkctmp(*) 
      integer , intent(in) :: iphead(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double)  :: wate(itdim,*) 
      real(double) , intent(in) :: bxv(*) 
      real(double) , intent(in) :: byv(*) 
      real(double) , intent(in) :: bzv(*) 
      real(double) , intent(inout) :: px(0:npart) 
      real(double) , intent(inout) :: py(0:npart) 
      real(double) , intent(inout) :: pz(0:npart) 
      real(double) , intent(inout) :: up(0:npart) 
      real(double) , intent(inout) :: vp(0:npart) 
      real(double) , intent(inout) :: wp(0:npart) 
      real(double) , intent(in) :: pxi(0:npart) 
      real(double) , intent(in) :: peta(0:npart) 
      real(double) , intent(in) :: pzta(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk_host, np, inew, jnew, knew, ijk_chek 
      real(double) :: normx, normy, normz, tsi, eta, dx1, dx2, dy1, dy2, dz1, &
         dz2, x_intersect, y_intersect, z_intersect, bxp, byp, bzp, rnorm, &
         dxdotn, udotb, udotn 
!-----------------------------------------------
!
!
!
!dir$ ivdep
!
      do n = 1, ncellsp 
!
         ijk_host = ijkctmp(n) 
         np = iphead(ijk_host) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
         ijk_chek = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
         tsi = pxi(np) - inew 
         eta = peta(np) - jnew 
!
         dx1 = (1. - eta)*(x(ijk_chek+iwid+kwid)-x(ijk_chek+kwid)) + eta*(x(&
            ijk_chek+iwid+jwid+kwid)-x(ijk_chek+jwid+kwid)) 
!
         dx2 = (1. - tsi)*(x(ijk_chek+jwid+kwid)-x(ijk_chek+kwid)) + tsi*(x(&
            ijk_chek+iwid+jwid+kwid)-x(ijk_chek+iwid+kwid)) 
!
         dy1 = (1. - eta)*(y(ijk_chek+iwid+kwid)-y(ijk_chek+kwid)) + eta*(y(&
            ijk_chek+iwid+jwid+kwid)-y(ijk_chek+jwid+kwid)) 
!
         dy2 = (1. - tsi)*(y(ijk_chek+jwid+kwid)-y(ijk_chek+kwid)) + tsi*(y(&
            ijk_chek+iwid+jwid+kwid)-y(ijk_chek+iwid+kwid)) 
!
         dz1 = (1. - eta)*(z(ijk_chek+iwid+kwid)-z(ijk_chek+kwid)) + eta*(z(&
            ijk_chek+iwid+jwid+kwid)-z(ijk_chek+jwid+kwid)) 
!
         dz2 = (1. - tsi)*(z(ijk_chek+jwid+kwid)-z(ijk_chek+kwid)) + tsi*(z(&
            ijk_chek+iwid+jwid+kwid)-z(ijk_chek+iwid+kwid)) 
!
         normx = dy1*dz2 - dy2*dz1 
         normy = dz1*dx2 - dz2*dx1 
         normz = dx1*dy2 - dx2*dy1 
!
         x_intersect = tsi*(x(ijk_chek+iwid+kwid)*(1.-eta)+x(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(x(ijk_chek+kwid)*(1.-eta)+x(ijk_chek+jwid&
            +kwid)*eta) 
!
         y_intersect = tsi*(y(ijk_chek+iwid+kwid)*(1.-eta)+y(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(y(ijk_chek+kwid)*(1.-eta)+y(ijk_chek+jwid&
            +kwid)*eta) 
!
         z_intersect = tsi*(z(ijk_chek+iwid+kwid)*(1.-eta)+z(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(z(ijk_chek+kwid)*(1.-eta)+z(ijk_chek+jwid&
            +kwid)*eta) 
!
         bxp = tsi*(bxv(ijk_chek+iwid+kwid)*(1.-eta)+bxv(ijk_chek+iwid+jwid+&
            kwid)*eta) + (1. - tsi)*(bxv(ijk_chek+kwid)*(1.-eta)+bxv(ijk_chek+&
            jwid+kwid)*eta) 
!
         byp = tsi*(byv(ijk_chek+iwid+kwid)*(1.-eta)+byv(ijk_chek+iwid+jwid+&
            kwid)*eta) + (1. - tsi)*(byv(ijk_chek+kwid)*(1.-eta)+byv(ijk_chek+&
            jwid+kwid)*eta) 
!
         bzp = tsi*(bzv(ijk_chek+iwid+kwid)*(1.-eta)+bzv(ijk_chek+iwid+jwid+&
            kwid)*eta) + (1. - tsi)*(bzv(ijk_chek+kwid)*(1.-eta)+bzv(ijk_chek+&
            jwid+kwid)*eta) 
!
!
!
         rnorm = 1./(normx**2 + normy**2 + normz**2) 
!
!
         dxdotn = ((px(np)-x_intersect)*normx+(py(np)-y_intersect)*normy+(pz(np&
            )-z_intersect)*normz)*rnorm 
!
         px(np) = px(np) - 2.*normx*dxdotn 
         py(np) = py(np) - 2.*normy*dxdotn 
         pz(np) = pz(np) - 2.*normz*dxdotn 
!
         udotb = (up(np)*bxp+vp(np)*byp+wp(np)*bzp)/(bxp**2 + byp**2 + bzp**2&
             + 1.E-10) 
!
!par      if(udotb.ne.0.0) then
!par      up(np)=-up(np)+2.*bxp*udotb
!par      vp(np)=-vp(np)+2.*byp*udotb
!par      wp(np)=-wp(np)+2.*bzp*udotb
!parc
!par      else
!
         udotn = (up(np)*normx+vp(np)*normy+wp(np)*normz)*rnorm 
!
         up(np) = up(np) - 2.*udotn*normx 
         vp(np) = vp(np) - 2.*udotn*normy 
         wp(np) = wp(np) - 2.*udotn*normz 
!
!par      endif
!
      end do 
!
      return  
      end subroutine parrefl 
